# simplified version of StrategyUnitTheme::su_theme_pal()
#   see: https://github.com/The-Strategy-Unit/StrategyUnitTheme
su_theme_pal <- function(palette = "main", ...) {
  pal <- c(orange = "#f9bf07",
           red    = "#ec6555",
           blue   = "#5881c1")
  grDevices::colorRampPalette(pal, ...)
}

If you haven’t already got these packages installed you will need the following to complete this workshop:

install.packages(c("tidyverse", "simmer", "simmer.plot", "gridExtra", "TruncatedNormal"))

1 Building the model - part 1

1.1 Libraries needed

suppressPackageStartupMessages({
  library(tidyverse)
  library(simmer)
  library(simmer.plot)
  library(TruncatedNormal)
  library(writexl)
  library(gridExtra)
})

1.2 Create distribution functions

For the patients interarrival rate, the time taken for asking the screening questions and the time for the receptionists to make the decision on what appointment to book.

You can use helper functions

# generic truncated normal distributions looks like:
rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
## [1] 0.6213012
dist_helper <- function(rfn, ...) {
  function() rfn(1, ...)
}
dist_patient_interarrival <- dist_helper(rexp, 50 / 60)
dist_screening            <- dist_helper(rtnorm, 4, 1,    0, Inf)
dist_decision             <- dist_helper(rtnorm, 1, 0.25, 0, Inf)
dist_booking              <- dist_helper(rtnorm, 1, 0.25, 0, Inf)

Alternatively:

dist_patient_interarrival <- function() rexp(1, 50 / 60)
dist_screening            <- function() rtnorm(n = 1, mu = 4, sd = 1,    lb = 0, ub = Inf)
dist_decision             <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
dist_booking              <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
dist_contact_mode         <- function() sample(1:2, 1, FALSE, c(0.15, 0.85))
dist_booking_type         <- function() sample(1:3, 1, FALSE, c(0.3, 0.6, 0.1))

1.3 Creating the trajectories for the branches

Think about the flow diagram of the process to create the trajectories for the model, let’s break it down into manageable sections.

1.3.1 ECR branch

Create a trajectory for this branch called branch_ecr. Within the trajectory use set_attribute via_ecr to label the patient with a value of 1, so we will know what route they have come to analyse later. Then seize a receptionist ready for the next stage of the process.

branch_ecr <- trajectory("branch ecr") %>%
  set_attribute("via_ecr", 1) %>%
  seize("receptionist")
branch_ecr
## trajectory: branch ecr, 2 activities
## { Activity: SetAttribute | keys: [via_ecr], values: [1], global: 0, mod: N, init: 0 }
## { Activity: Seize        | resource: receptionist, amount: 1 }

1.3.2 Phone branch

create a trajectory for this branch called branch_phone. Within the trajectory use set_attribute via_ecr to label the patient with a value of 0 (as opposed to 1 which are ECR route patients).

The patients renege after 5 minutes so log this and set the attribute "reneged" to 1. Then seize a receptionist ready for the next stage. Once the patient has a receptionist abort the renege and add timeout of distribution screening

renege_hungup <- trajectory("hung up") %>%
  set_attribute("reneged", 1) %>%
  log_("lost my patience. Hanging up...")

branch_phone <- trajectory("branch phone") %>%
  set_attribute("via_ecr", 0) %>%
  renege_in(t = 5, out = renege_hungup) %>%
  seize("receptionist") %>%
  renege_abort() %>% # now have a receptionist the patient won't renege
  timeout(dist_screening)

1.3.3 Booking option branches

Receptionist is deciding if a patient needs phone appointment, needs phone appointment or does not need any appointment.

branch_booking_none <- trajectory("no appointment needed") %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 0)
branch_booking_none
## trajectory: no appointment needed, 2 activities
## { Activity: Release      | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [0], global: 0, mod: N, init: 0 }
branch_booking_phone <- trajectory("booking a phone") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 1)
branch_booking_phone
## trajectory: booking a phone, 3 activities
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [1], global: 0, mod: N, init: 0 }
branch_booking_f2f <- trajectory("booking a f2f") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 2)
branch_booking_f2f
## trajectory: booking a f2f, 3 activities
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [2], global: 0, mod: N, init: 0 }

1.4 Patient Trajectory

patient_flow <- trajectory("patient flow") %>%
  branch(dist_contact_mode, TRUE,
         branch_ecr,
         branch_phone) %>%
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none,
         branch_booking_phone,
         branch_booking_f2f)
plot(patient_flow, fill = su_theme_pal("main"))
# you can change the verbose argument to TRUE to get more information
plot(patient_flow, fill = su_theme_pal("main"), verbose = TRUE)

1.5 Add resources and generator to the environment

You can use to() function in the generator, so that it completes the arrivals. Use to() with caution as the stop time is that of the generator and not the simulation. The simulation will run until all arrivals are finished, so we are assuming that staff stay overtime

env_part1 <- simmer("part 1") %>%
  add_resource("receptionist", 2) %>%
  add_generator("patient", patient_flow, to(8 * 60, dist_patient_interarrival), mon = 2)

1.6 Run the model

1.7 Glimpse results to check run

env_part1 %>%
  get_mon_arrivals() %>%
  arrange(start_time) %>%
  head(10)
name start_time end_time activity_time finished replication
patient0 0.9062182 8.666809 7.760591 TRUE 1
patient1 1.0810663 2.880660 1.799594 TRUE 1
patient2 1.6043486 7.786732 4.906072 TRUE 1
patient3 2.0532035 7.053204 0.000000 TRUE 1
patient4 2.3693078 7.369308 0.000000 TRUE 1
patient5 2.7333249 7.733325 0.000000 TRUE 1
patient6 3.6477607 13.937513 6.150780 TRUE 1
patient7 3.9308550 13.606157 4.939348 TRUE 1
patient8 4.9175231 9.917523 0.000000 TRUE 1
patient9 7.7549414 12.754941 0.000000 TRUE 1
env_part1 %>%
  get_mon_arrivals() %>%
  arrange(start_time) %>%
  tail(10)
name start_time end_time activity_time finished replication
389 patient388 472.4713 477.9166 5.445350 TRUE 1
390 patient389 472.8399 483.2193 6.273473 TRUE 1
391 patient390 473.1776 481.5829 3.666322 TRUE 1
392 patient391 474.4332 479.4332 0.000000 TRUE 1
393 patient392 475.5632 480.5632 0.000000 TRUE 1
394 patient393 476.6170 487.7046 6.121711 TRUE 1
395 patient394 476.7175 481.7175 0.000000 TRUE 1
396 patient395 477.5498 482.5498 0.000000 TRUE 1
397 patient396 478.8235 490.1992 6.979876 TRUE 1
398 patient397 479.2955 484.2955 0.000000 TRUE 1

1.8 Questions

If you are new to R and struggling with data manipulation, you can export the results as xls or csv. Or email us and we can help you!

write_xlsx(get_mon_arrivals(env_part1), "arrivals_part1.xlsx")
write_xlsx(get_mon_resources(env_part1), "resources_part1.xlsx")
write_xlsx(get_mon_attributes(env_part1), "attributes_part1.xlsx")

1.8.1 Count reneged

env_part1 %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  summarise(countreneged = n())
countreneged
179

n = 148 with 2 receptionists

1.8.2 Count F2F

env_part1 %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  summarise(countF2F = n())
countF2F
30

n = 27 with 2 receptionists

1.8.3 Call waits

Get attributes of finished patients from get_mon_attributes(), join with get_mon_arrivals(). Because of varying number of decimal places we get some rounding errors resulting in near zero but negative results. That’s why when we are doing data wrangling, wee use pmax (), so it takes 0 or the wait time which ever is greater. We can then calculate summary statistics.

part1_attributes <- env_part1 %>%
  get_mon_attributes() %>%
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)

part1_callwait_results <- env_part1 %>%
  get_mon_arrivals(ongoing = TRUE) %>% # not needed as there aren't any due to the generator and run()
  inner_join(part1_attributes, by = c("name", "replication")) %>%
  filter(via_ecr == 0) %>%
  mutate(wait_time_call = case_when(
    # the caller reneged
    reneged == 1 ~ 5,
    # in queue still - these could also still be in progress
    # (would need to know time resource seized to calculate wait time),
    # except we used to() in generator and run() so we have no unfinished
    # !finished ~ (10 * 60) - start_time,
    # everything else, (assume has been answered)
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

part1_callwait_summary <- part1_callwait_results %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished))

part1_callwait_summary
median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished
5 4.223388 1.733617 0 5 3.906749 5 0 329

with 2 receptionists

#     median_wait   mean_wait      var_wait    min_wait    max_wait      Qtile25th_wait    Qtile75th_wait
#1    4.839153        3.811972      2.657451         0         5           3.010095              5

#     countunfinished countfinished
#1               0           334

2 Scenario additional model features

Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog - let’s look at 1 replication

attributes <- env_part1 %>%
  get_mon_attributes() %>%
  filter(replication == 1)

test <- env_part1 %>%
  get_mon_arrivals() %>%
  filter(finished == FALSE, replication == 1) %>%
  inner_join(attributes, by = "name") %>%
  arrange(start_time)

test
name start_time end_time activity_time finished replication.x time key value replication.y
env_part1 %>%
  get_mon_attributes() %>%
  filter(replication == 1, key == "via_ecr", value == 1) %>%
  count()
n
69

You can see that everyone finished and we had total of 66 ECR.

3 Extra model features - part 2

3.1 New environment

Create a new environment

env_part2 <- simmer("part2")
env_part2
## simmer environment: part2 | now: 0 | next: 
## { Monitor: in memory }

3.2 Add time-dependent interarrival times

To include time dependent interarrivals create a different distribution, one for each time period (3 in total). You may also find using these helper functions to switch between minutes, hours and days useful:

t_minute <- function(t) t
t_hour   <- function(t) t_minute(t) * 60
t_day    <- function(t) t_hour(t) * 24

dist_patient_interarrival1 <- function() rexp(1, 120 / t_hour(2))
dist_patient_interarrival2 <- function() rexp(1, 240 / t_hour(6))
dist_patient_interarrival3 <- function() rexp(1,  40 / t_hour(2))

3.3 Phone branch with prioritisation

We are also going prioritise the patients contacting the surgery via the phone.

In the ECR branch, no changes are needed. The default priority for all arrivals is set up to 0.

In the phone branch, we want to set a higher priority using set_prioritization(), that receptionists prioritise answering calls over reviewing electronic requests but without interrupting current tasks (prioirity = 1, preemptive = 1, restart = FALSE). If you didn’t use now(env_part1) to set attributes you can use join(branch_phone) instead to link after the prioritisation back to your original phone branch.

branch_phone2 <- trajectory("branch phone prioritised") %>%
  set_prioritization(c(1, 1, FALSE)) %>%
  join(branch_phone)

plot(branch_phone2, fill = su_theme_pal("main"))

from the help page for add_generator():

  • priority: the priority of each arrival (a higher integer equals higher priority; defaults to the minimum priority, which is 0.

  • preemptible: if a seize occurs in a preemptive resource, this parameter establishes the minimum incoming priority that can preempt these arrivals (an arrival with a priority greater than preemptible gains the resource). In any case, preemptible must be equal or greater than priority, and thus only higher priority arrivals can trigger preemption.

  • restart: whether the activity must be restarted after being preempted.

3.4 Patient Trajectory

same as Part 1 but with the new prioritised phone branch

patient_flow2 <- trajectory("patient flow2") %>%
  branch(dist_contact_mode, TRUE,
         branch_ecr,
         branch_phone2) %>% #<-- changed the phone branch
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none,
         branch_booking_phone,
         branch_booking_f2f)
plot(patient_flow2, fill = su_theme_pal())

3.5 Schedule for staffing

Use the schedule(timetable, values, period) function and set it to repeat the schedule every 24 hour period. Work out from the numbers on each shift what times capacity changes and how many are working at that point:

e.g. 0-8am = 0, 8-9am = 1, 9-9.30 am etc. becomes:

sch_receptionists <- schedule(
  timetable = t_hour(c(0.0, 8.0, 9.0, 9.5, 12.0, 15.5, 18.0)),
  values    =        c(0,   1,   2,   4,    3,    1,    0),
  period    = t_hour(24)
)

3.6 Add resources and generator to the environment

Add the receptionists as a resource using the schedule you created for the capacity. Either add three different generators to the model environment setting the times they operate with from_to or use functional programming to control this. You can look back at the training slides and examples to help you.

env_part2 %>%
  add_resource("receptionist", capacity = sch_receptionists) %>%
  add_generator("patient[8am-10am]", patient_flow2,
                from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[10am-4pm]", patient_flow2,
                from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[4pm-6pm]", patient_flow2,
                from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
                mon = 2)
## simmer environment: part2 | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }

3.7 Run model

Set the random seed to 1 to give us the same results each time. Run your model for 5 days (52460) or use one of the time helper functions)

3.8 Questions

Repeat your analysis of call waits, reneging and F2F appointments as you did for part 1 using the environment for part 2. Plot utilisation of the resources using the plot function with get_mon_arrivals() and the metric “waiting_time”.

3.8.1 Count reneged

env_part2 %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  summarise(countreneged = n())
countreneged
582

n = 582 with 2 receptionists

3.8.2 Count F2F

env_part2 %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  summarise(countF2F = n())
countF2F
160

n = 160 with 2 receptionists

3.8.3 Call waits

part2_attributes <- env_part2 %>%
  get_mon_attributes() %>%
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)

part2_callwait_results <- env_part2 %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  full_join(part2_attributes, by = c("name", "replication")) %>% # or can use inner_join
  filter(via_ecr == 0) %>%
  mutate(wait_time_call = case_when(
    reneged == 1 ~ 5,
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

part2_callwait_summary <- part2_callwait_results %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished))

part2_callwait_summary
median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished
3.825883 3.183036 3.479728 0 5 1.464796 5 0 1726
#   median_wait     mean_wait     var_wait    min_wait    max_wait    Qtile25th_wait    Qtile75th_wait
#1    3.825883      3.183036      3.479728          0          5          1.464796              5
#   countunfinished countfinished
#1               0          1726

3.9 What else would you do to improve this model?

Some ideas:

  • prevent staff staying to finish all just the one they are on,
  • multiple runs,
  • staff breaks, availability, are they also doing other tasks?

3.10 What other metrics would you want to measure / the practice find useful?

Below are some ideas on additional analysis you might want to do

Plot waiting time over time

plot(get_mon_arrivals(env_part2), metric = "waiting_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

boxplot to show variation in metrics e.g. call_wait

g1 <- ggplot() +
  geom_boxplot(data = part2_callwait_results,
               aes(x = "rep=1", y = wait_time_call),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  geom_point(data = part2_callwait_summary,
             aes(x = "rep=1", y = mean_wait),
             shape = 19,
             size = 2,
             colour = "red") +
  labs(x = "Part2",
       y = "call wait (mins)",
       title = "Boxplot of call wait from single replication") +
  theme(axis.text = element_text(size = 11))

g1

Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog:

attributes <- env_part2 %>%
  get_mon_attributes() %>%
  filter(replication == 1)

test <- env_part2 %>%
  get_mon_arrivals() %>%
  filter(finished == FALSE, replication == 1) %>%
  inner_join(attributes, by = "name") %>%
  arrange(start_time)

test %>%
  filter(key == "via_ecr", value != 1)
name start_time end_time activity_time finished replication.x time key value replication.y
env_part2 %>%
  get_mon_attributes() %>%
  filter(replication == 1, key == "via_ecr", value == 1) %>%
  count()
n
357
test %>%
  head(10) %>%
  as_tibble()
name start_time end_time activity_time finished replication.x time key value replication.y
patient[4pm-6pm]191 6774.405 NA NA FALSE 1 6774.405 via_ecr 1 1
patient[4pm-6pm]192 6775.585 NA NA FALSE 1 6775.585 via_ecr 1 1
patient[4pm-6pm]201 6790.286 NA NA FALSE 1 6790.286 via_ecr 1 1

All unfinished are ECR, because calls were prioritised ahead of them. 357 total ECR.

# A tibble: 3 x 10 - 3 are unfinished in a week out of 357
# name                start_time end_time activity_time finished replication.x  time key     value replication.y
#<chr>                    <dbl>    <dbl>         <dbl> <lgl>            <int> <dbl> <chr>   <dbl>         <int>
#  patient[4pm-6pm]191      6774.       NA            NA FALSE                1 6774. via_ecr     1             1
#  patient[4pm-6pm]192      6776.       NA            NA FALSE                1 6776. via_ecr     1             1
#  patient[4pm-6pm]201      6790.       NA            NA FALSE                1 6790. via_ecr     1             1

Look at resources (reception staff in this example)

resources2 <- get_mon_resources(env_part2)
plot(resources2, metric = "usage", "receptionist", items = "server", steps = TRUE)

plot(resources2, metric = "utilization", "receptionist")

# show the top 10 results
head(resources2, n = 10)
resource time server queue capacity queue_size system limit replication
receptionist 480.0000 0 0 1 Inf 0 Inf 1
receptionist 480.0000 1 0 1 Inf 1 Inf 1
receptionist 481.1816 1 1 1 Inf 2 Inf 1
receptionist 481.6397 1 2 1 Inf 3 Inf 1
receptionist 484.5346 1 3 1 Inf 4 Inf 1
receptionist 484.9699 1 4 1 Inf 5 Inf 1
receptionist 485.1100 1 3 1 Inf 4 Inf 1
receptionist 486.1832 1 4 1 Inf 5 Inf 1
receptionist 486.6397 1 3 1 Inf 4 Inf 1
receptionist 487.5739 1 4 1 Inf 5 Inf 1
# list of hours
HrList <- data.frame(
  resource = "receptionist",
  time = seq(from = 0, to = t_day(5), by = 60),
  server = NA,
  queue = NA,
  capacity = NA,
  queue_size = Inf,
  system = NA,
  limit = Inf,
  replication = 1
)

# which aren't in the data
AddSnaps <- filter(HrList, !(time %in% resources2$time))

resources2_hourly <- bind_rows(resources2, AddSnaps) %>%
  arrange(time) %>%
  mutate(server = ifelse(is.na(server), lag(server), server),
         queue = ifelse(is.na(queue), lag(queue), queue),
         capacity = ifelse(is.na(capacity), lag(capacity), capacity),
         system = ifelse(is.na(system), lag(system), system))

resourceutilisation_byhour <- resources2_hourly %>%
  arrange(time) %>%
  mutate(slot_hour = floor(time / 60)) %>%
  mutate(slot_day = ceiling(slot_hour / 24)) %>%
  mutate(slot = slot_hour - ((slot_day - 1) * 24)) %>%
  group_by(slot) %>%
  mutate(dt = time - lag(time),
         capacity = ifelse(capacity < server, server, capacity),
         in_use = dt * lag(server / capacity)) %>%
  summarise(utilised_percent = sum(in_use, na.rm = TRUE) / sum(dt, na.rm = TRUE) * 100)

resourceutilisation_byhour
slot utilised_percent
1 0.0000000
2 0.0000000
3 0.0000000
4 0.0000000
5 0.0000000
6 0.0000000
7 0.0000000
8 100.0000000
9 100.0000000
10 99.9259080
11 63.7850307
12 99.6780740
13 91.7546945
14 75.7575154
15 99.8493045
16 99.8748793
17 100.0000000
18 0.1985734
19 0.0000000
20 0.0000000
21 0.0000000
22 0.0000000
23 0.0000000
24 0.0000000

3.11 Alternative way to calculate waiting time - option B

Option B. If you had used now(env_part1) to set attributes then you will need to copy your code and change the environment for all parts of the trajectory in which it was used. You then need to create the overall trajectory again using the new version of the phone branch (and other parts if you have used now()).

Firstly, create a new environment.

env_part2B <- simmer("part 2B")
env_part2B
## simmer environment: part 2B | now: 0 | next: 
## { Monitor: in memory }

Change references to part env_part2B in all trajectories

renege_hungup2B <- trajectory("hung up") %>%
  set_attribute("reneged", 1) %>%
  set_attribute("reneged_at", function() now(env_part2B)) %>% # changed to env_part2B
  log_("lost my patience. Hanging up...")

branch_phone2B <- trajectory("branch phone") %>%
  set_prioritization(c(1, 1, FALSE)) %>% # set to increase priority from default of 0 to 1, no preempt and no restart
  set_attribute("via_ecr", 0) %>%
  set_attribute("arrival", function() now(env_part2B)) %>% # changed to env_part2B
  renege_in(t = 5, out = renege_hungup) %>%
  seize("receptionist") %>%
  set_attribute("call_taken", function() now(env_part2B)) %>% # changed to env_part2B
  renege_abort() %>%
  timeout(dist_screening)

branch_booking_none2B <- trajectory("no appointment needed") %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 0) %>%
  set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B

branch_booking_phone2B <- trajectory("booking a f2f") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 1) %>%
  set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B

branch_booking_f2f2B <- trajectory("booking a f2f") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 2) %>%
  set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B

patient_flow2B <- trajectory("patient flow2B") %>%
  branch(dist_contact_mode, TRUE,
         branch_ecr,
         branch_phone2B) %>%
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none2B,
         branch_booking_phone2B,
         branch_booking_f2f2B)

Populate environment using the trajectory for patient flow 2B

env_part2B %>%
  add_resource("receptionist", capacity = sch_receptionists) %>%
  add_generator("patient[8am-10am]", patient_flow2B,
                from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[10am-4pm]", patient_flow2B,
                from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[4pm-6pm]", patient_flow2B,
                from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
                mon = 2)
## simmer environment: part 2B | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }

Run your model for 5 days

3.11.1 Part 2 Questions for Option B

What is the average wait for calls to be answered, number abandoning their calls (hanging up) and the number booked for a face to face appointment in this simulated week?

env_part2B %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  summarise(countreneged = n())
countreneged
582

n=582

How many people get a F2F booking?

env_part2B %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  summarise(countF2F = n())
countF2F
160

n=160

What are the call waits?

part2B_callwait_results <- env_part2B %>%
  get_mon_attributes() %>%
  pivot_wider(-c(time, replication), names_from = key, values_from = value) %>%
  filter(via_ecr == 0, arrival >= 0) %>%
  mutate(wait_time_call = case_when(
    !is.na(reneged) ~ 5,
    is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
    TRUE ~ call_taken - arrival # has been answered
  ))

part2B_callwait_summary <- part2B_callwait_results %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(is.na(call_end) & is.na(reneged)),
            countfinished = sum(!is.na(call_end) | !is.na(reneged)))

part2B_callwait_summary
median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished
3.825883 3.183036 3.479728 0 5 1.464796 5 0 1726
#   # A tibble: 1 x 9
#       median_wait     mean_wait     var_wait  min_wait    max_wait Qtile25th_wait Qtile75th_wait
#1         3.83          3.18        3.48        0              5           1.46              5
# countunfinished countfinished
# 0                 1726

4 Multiple runs and Scenarios - part 3

Now we are going to run our model 50 times and analyse different scenarios.

4.1 Replication of baseline

Use lapply to loop through and replicate the environment you created in Part 2 100 times

4.2 get baseline results and plots

Update some of the code you used to look at results in pasrt 1 and 2. This time group_by(replication) before you summarise. E.g count of the reneged, then ungroup(). This will give the result for each replication. It is these that are then summarised again to give your final results for i.e. the mean of the mean number reneged in each replication.

4.2.1 Count reneged

envs_baseline %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))
countreneged_mean
579.86

n=580

4.2.2 Count F2F

envs_baseline %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))
countF2F_mean
143.04

n=143

Before we look at call waits let’s understand a little more about unfinished. using the staff schedule means that the final staff member finishes at 6pm on the 5th day so at (42460)+(18*60)= 6840, there are 4 in the system, 1 patient in progress and 3 waiting the staff member finishes that one at 6845 and then leaves with 3 remaining unfinished a similar thing occurs at the end of day 3:

get_mon_resources(envs_baseline) %>%
  filter(replication == 1) %>%
  tail(n = 10) # show last 10 rows
resource time server queue capacity queue_size system limit replication
4184 receptionist 6824.155 1 4 1 Inf 5 Inf 1
4185 receptionist 6827.207 1 5 1 Inf 6 Inf 1
4186 receptionist 6829.155 1 4 1 Inf 5 Inf 1
4187 receptionist 6829.959 1 3 1 Inf 4 Inf 1
4188 receptionist 6832.803 1 4 1 Inf 5 Inf 1
4189 receptionist 6835.244 1 3 1 Inf 4 Inf 1
4190 receptionist 6835.407 1 4 1 Inf 5 Inf 1
4191 receptionist 6839.038 1 3 1 Inf 4 Inf 1
4192 receptionist 6840.000 1 3 0 Inf 4 Inf 1
4193 receptionist 6845.235 0 3 0 Inf 3 Inf 1
get_mon_resources(envs_baseline) %>%
  filter(replication == 1, time < 6240) %>%
  tail(n = 10)
resource time server queue capacity queue_size system limit replication
3326 receptionist 5387.154 1 4 1 Inf 5 Inf 1
3327 receptionist 5387.827 1 3 1 Inf 4 Inf 1
3328 receptionist 5391.761 1 4 1 Inf 5 Inf 1
3329 receptionist 5392.154 1 3 1 Inf 4 Inf 1
3330 receptionist 5392.374 1 4 1 Inf 5 Inf 1
3331 receptionist 5392.431 1 5 1 Inf 6 Inf 1
3332 receptionist 5393.999 1 4 1 Inf 5 Inf 1
3333 receptionist 5397.431 1 3 1 Inf 4 Inf 1
3334 receptionist 5400.000 1 3 0 Inf 4 Inf 1
3335 receptionist 5401.615 0 3 0 Inf 3 Inf 1
envs_baseline %>%
  get_mon_arrivals(ongoing = TRUE) %>% # 101972 rows
  filter(replication == 1, is.na(end_time))
name start_time end_time activity_time finished replication
patient[4pm-6pm]212 -1.000 NA NA FALSE 1
patient[4pm-6pm]201 6790.286 NA NA FALSE 1
patient[4pm-6pm]191 6774.405 NA NA FALSE 1
patient[4pm-6pm]192 6775.585 NA NA FALSE 1
patient[10am-4pm]1218 -1.000 NA NA FALSE 1
patient[8am-10am]653 -1.000 NA NA FALSE 1

4.2.3 Call waits

baseline_attributes <- envs_baseline %>%
  get_mon_attributes() %>% #202732 obs
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
#101468 obs

baseline_callwait_results <- envs_baseline %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  distinct() %>% # goes down to 101618, removes duplicates of unfinished
  # still more in arrivals than attributes as attributes only has finished
  full_join(baseline_attributes, by = c("name", "replication")) %>%
  filter(via_ecr == 0) %>% # phone only also removes those yet to start e.g. start_time -1
  mutate(wait_time_call = case_when(
    # the caller reneged
    reneged == 1 ~ 5,
    # unfinished still in the system,
    !finished ~ t_day(5) - start_time,
    # everything else, (assume has been answered)
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

## calculate summary statistics
baseline_callwait_replication_summary <- baseline_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished)) %>%
  mutate(scenario = "baseline")

# show top 10 rows
head(baseline_callwait_replication_summary, n = 10)
replication median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished scenario
1 3.825883 3.183036 3.479728 0 5 1.464796 5 0 1726 baseline
2 3.744117 3.221769 3.258414 0 5 1.555424 5 0 1770 baseline
3 3.775021 3.116034 3.682359 0 5 1.206570 5 0 1772 baseline
4 3.952666 3.235613 3.446565 0 5 1.535425 5 0 1722 baseline
5 3.810451 3.201863 3.433046 0 5 1.446567 5 0 1736 baseline
6 3.550199 3.053321 3.688564 0 5 1.141632 5 0 1699 baseline
7 3.525437 3.005393 3.874786 0 5 1.005038 5 0 1652 baseline
8 3.808467 3.141373 3.616938 0 5 1.237069 5 0 1733 baseline
9 3.892609 3.168985 3.635974 0 5 1.225592 5 0 1750 baseline
10 3.623343 3.068928 3.714392 0 5 1.107349 5 0 1662 baseline
baseline_callwait_summary <- baseline_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "baseline") # to identify it later

baseline_callwait_summary
median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished scenario
3.765179 3.147497 3.580269 0 5 1.334405 5 0 1722.8 baseline

5 Scenario

If a neighbouring practice closed leading to an increase in demand of 45%. Whilst the practice were able to better advertise the ECR route to patients such that 40% of all the patients (old and new) were able to use this route. What impact would this have on the call wait and reneging of patients?

5.1 create patient flow 3 for this scenario with 40% ECR

dist_contact_mode_S <- function() sample(1:2, 1, FALSE, c(0.4, 0.6))
patient_flow3 <- trajectory("patient flow3") %>%
  branch(dist_contact_mode_S, TRUE, # <- changed to the scenario distribution
         branch_ecr,
         branch_phone2) %>%
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none,
         branch_booking_phone,
         branch_booking_f2f)

5.2 scenario interarrivals

dist_patient_interarrival1_S <- function() rexp(1, (120 * 1.45) / t_hour(2))
dist_patient_interarrival2_S <- function() rexp(1, (240 * 1.45) / t_hour(6))
dist_patient_interarrival3_S <- function() rexp(1,  (40 * 1.45) / t_hour(2))

5.3 lapply as before but using patient_flow3 and increasing frquency of arrivals

5.4 Scenario results

Repurpose some of the code you used to look at results in parts 1 and 2 of the workshop. Remember to store the results with a different name to your baseline results.

5.4.1 Count reneged

envs_scenario %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))
countreneged_mean
603.16

n=752

5.4.2 Count F2F

envs_scenario %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))
countF2F_mean
171.56

n=143

5.4.3 Call waits

scenario_attributes <- envs_scenario %>%
  get_mon_attributes() %>%
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)

scenario_callwait_results <- envs_scenario %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  distinct() %>%
  full_join(scenario_attributes, by = c("name", "replication")) %>%
  filter(via_ecr == 0) %>%
  mutate(wait_time_call = case_when(
    # the caller reneged
    reneged == 1 ~ 5,
    # unfinished still in the system,
    !finished ~ t_day(5) - start_time,
    # everything else, (assume has been answered)
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

## calculate summary statistics
scenario_callwait_replication_summary <- scenario_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished)) %>%
  mutate(scenario = "scenario")

scenario_callwait_summary <- scenario_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "scenario")

# show top 10 results
head(scenario_callwait_summary, n = 10)
median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished scenario
3.890618 3.289315 3.122619 0.001257 5 1.598594 5 0 1747.56 scenario

5.5 T-test

Use t.test() function to compare waiting times between baseline and scenario and see if there is a significant difference in the results.

5.5.1 Version taught in slidepack:

t.test(scenario_callwait_results$wait_time_call, baseline_callwait_results$wait_time_call)
## 
##  Welch Two Sample t-test
## 
## data:  scenario_callwait_results$wait_time_call and baseline_callwait_results$wait_time_call
## t = 16.071, df = 172339, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1241180 0.1585979
## sample estimates:
## mean of x mean of y 
##  3.290622  3.149264

5.5.2 Paired T-test

This could actually be a paired t-test because they use the same common random number set and have the same number of replications therefore replication 1 of baseline and scenario are correlated (unless the scenario itself changes the andom number sampling). we should check the variances of the differences is less than the individual variances before using paired t test. To do it, we need to:

  • combine a dataset
  • compute the difference between the scenario and baseline mean waits for each pair of replications
  • conduct Shapiro-Wilk normality test for the differences and check p-value
  • If p-value is >0.05, we conduct paired t-test
part3_summary_byreplication <- bind_rows(
  baseline_callwait_replication_summary,
  scenario_callwait_replication_summary
)

d <- with(
  part3_summary_byreplication,
  mean_wait[scenario == "baseline"] - mean_wait[scenario == "scenario"]
)
shapiro.test(d) # => p-value = 0.6375
## 
##  Shapiro-Wilk normality test
## 
## data:  d
## W = 0.95276, p-value = 0.04435
res <- t.test(mean_wait ~ scenario, data = part3_summary_byreplication, paired = TRUE)
res
## 
##  Paired t-test
## 
## data:  mean_wait by scenario
## t = -8.5021, df = 49, p-value = 3.295e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1753380 -0.1082975
## sample estimates:
## mean of the differences 
##              -0.1418177

p-value is greater than the significance level 0.05 implying that the distribution of the differences (d) are not significantly different from normal distribution. therefore we do a paired T-Test

data:  mean_wait by scenario
t = -26.092, df = 49, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  -0.4228384 -0.3623634
sample estimates:
  mean of the differences
-0.3926009

The p-value of the test is < 2.2e-16, which is less than the significance level alpha = 0.05. We can then reject null hypothesis and conclude that the average wait in the baseline is significantly different from the average wait in the scenario

5.5.3 Boxplot of the call waits in the baseline vs scenario

Using the ggplot2 package and geom_boxplot(), create boxplots of the means of the replications of the call waiting time in the baseline compared to the scenario. Set theme and title if you wish. To compare boxplots, you can either

  1. create 2 different plots and put them on one chart using gridExtra package
a <- ggplot(baseline_callwait_replication_summary) +
  geom_boxplot(aes(x = scenario, y = mean_wait),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  labs(y = "mean call wait (mins)") +
  ylim(0, 4) +
  theme(axis.title.x = element_blank(),
        axis.text = element_text(size = 11))

b <- ggplot(scenario_callwait_replication_summary) +
  geom_boxplot(aes(x = scenario, y = mean_wait),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  labs(y = "mean call wait (mins)") +
  ylim(0, 4) +
  theme(axis.title.x = element_blank(),
        axis.text = element_text(size = 11))

grid.arrange(a, b, ncol = 2, nrow = 1, top = "Call wait mean results of baseline(left) and scenario(right)")

  1. create one plot from combined dataset. Using bind_rows(), create one dataset with both scenarios and plot them together using ggplot2 already have the averages of by replications in part3_summary_byreplication
part3_summary <- bind_rows(baseline_callwait_summary, scenario_callwait_summary)

# plot
g <- ggplot() +
  geom_boxplot(data = part3_summary_byreplication,
               aes(x = scenario, y = mean_wait),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  geom_point(data = part3_summary,
             aes(x = scenario, y = mean_wait),
             shape = 19,
             size = 2,
             colour = "red") +
  labs(x = "scenario", y = "mean call wait (mins)", title = "Boxplot of mean call wait from 50 replications") +
  theme(axis.text = element_text(size = 11))

g

5.5.4 Part 3 Questions

Looking at the results of baseline of 50 replications, what can you say about the results?

What do you interpret from the results of the scenario in comparison to the baseline?

What other analysis might you want to do of the model results?

Additional analysis: example if you want to compare ECR to Phones waits might do it this way:

baseline_wait_results <- envs_baseline %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  distinct() %>%
  full_join(baseline_attributes, by = c("name", "replication")) %>%
  filter(start_time >= 0) %>% # so use this to remove the not started e.g. start_time -1
  mutate(wait_time_call = case_when(
    reneged == 1 ~ 5,
    !finished ~ t_day(5) - start_time,
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

baseline_wait_replication_summary <- baseline_wait_results %>%
  group_by(replication, via_ecr) %>% #group by both of these to get results for phone and ECR
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished)) %>%
  mutate(scenario = "baseline")
## `summarise()` has grouped output by 'replication'. You can override using the `.groups` argument.
# show top 10 results
head(baseline_wait_replication_summary, n = 10)
replication via_ecr median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished scenario
1 0 3.825883 3.183036 3.479728 0 5.000 1.464796 5.00000 0 1726 baseline
1 1 56.571437 121.432726 60680.123013 0 1129.364 14.716072 99.92010 3 354 baseline
2 0 3.744117 3.221769 3.258414 0 5.000 1.555424 5.00000 0 1770 baseline
2 1 65.218754 149.913667 69329.798636 0 1080.604 15.965540 149.12634 1 291 baseline
3 0 3.775021 3.116034 3.682359 0 5.000 1.206570 5.00000 0 1772 baseline
3 1 32.526169 64.839190 22473.421606 0 1075.412 5.002446 64.15030 6 314 baseline
4 0 3.952666 3.235613 3.446565 0 5.000 1.535425 5.00000 0 1722 baseline
4 1 41.584256 83.401127 35084.985697 0 1099.012 12.369562 73.55387 1 303 baseline
5 0 3.810451 3.201863 3.433046 0 5.000 1.446567 5.00000 0 1736 baseline
5 1 40.110651 76.963237 27211.750428 0 1008.028 9.223911 84.78485 0 299 baseline
baseline_wait_summary <- baseline_wait_replication_summary %>%
  group_by(via_ecr) %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished))

baseline_wait_summary
via_ecr median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished
0 3.765179 3.147497 3.580269 0 5.000 1.334405 5.00000 0.00 1722.80
1 42.015194 94.007905 41412.497885 0 1052.733 10.070932 85.24734 4.08 302.48

What other ways might you want to present data?

For example, plot results from all replications

plot(get_mon_arrivals(envs_baseline), metric = "waiting_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 1062 rows containing non-finite values (stat_smooth).
## Warning: Removed 1062 row(s) containing missing values (geom_path).

plot(get_mon_arrivals(envs_scenario), metric = "waiting_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 29002 rows containing non-finite values (stat_smooth).
## Warning: Removed 29002 row(s) containing missing values (geom_path).

Histogram

ggplot(part3_summary_byreplication, aes(mean_wait)) +
  geom_histogram(fill = "orange") +
  labs(title = "Distribution of average call wait in 50 replications") +
  facet_wrap(vars(scenario), ncol = 2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

5.5.5 Additional questions if you want to explore model further

What other opportunities are there for improvement of the process?

What other scenarios might you want to consider with the practice? What would your performance indicators be? How would you incorporate them in the model? Have a go if you have time. warm up period, staff rota, percentage coming via ecr, include preemption and restart

What else would need to be considered before they could implement changes suggested by the model? costs, other tasks the staff do, feasibility of suggested rota in real life

Are we running the model for a sufficient number of runs? process converted to R from: Simulations: The Practice of Model Development and Use, Robinson, S., 2nd Ed. 2014, Palgrave Macmillan

Remember back to our t test results we saved in res gave a confidence interval based on student’s t distribution

res$conf.int[1]
## [1] -0.175338
res$conf.int[2]
## [1] -0.1082975

Create a blank data frame to hold the cumulative results

check_baseline_replications <- data.frame(matrix(ncol = 5, nrow = 0))
x <- c("replication", "cmean", "tL", "tU", "percentdeviation")
colnames(check_baseline_replications) <- x

#for loop to see how the percent deviation of the confidence limits from the mean vary by the number of replications
for (i in 2:50) {
  df <- baseline_callwait_replication_summary %>%
    dplyr::select(replication, mean_wait) %>%
    filter(replication <= i) %>%
    summarise(replication = i,
              cmean = mean(mean_wait),
              tL = t.test(mean_wait)$conf.int[1],
              tU = t.test(mean_wait)$conf.int[2],
              percentdeviation = (cmean - tL) / cmean * 100) %>%
    data.frame()
  check_baseline_replications <- bind_rows(check_baseline_replications, df)
}

# show top 10 results
head(check_baseline_replications, n = 10)
replication cmean tL tU percentdeviation
2 3.202403 2.956327 3.448479 7.684104
3 3.173613 3.040728 3.306499 4.187198
4 3.189113 3.103887 3.274340 2.672417
5 3.191663 3.133636 3.249691 1.818099
6 3.168606 3.094869 3.242343 2.327118
7 3.145290 3.062986 3.227594 2.616748
8 3.144801 3.075910 3.213691 2.190611
9 3.147488 3.087915 3.207060 1.892702
10 3.139632 3.084423 3.194840 1.758443
11 3.153205 3.095464 3.210947 1.831200

6 Part 3 if you use now() function

To include the set_attribute using now() we reference the correct environment for each replication

We can now analyse results

6.0.1 Count reneged

envs_baselineB %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))
countreneged_mean
579.86

n=580

6.0.2 Count F2F

envs_baselineB %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))
countF2F_mean
143.04

n=143

6.0.3 Call waits

baselineB_callwait_results <- envs_baselineB %>%
  get_mon_attributes() %>%
  pivot_wider(-c(time), names_from = key, values_from = value) %>%
  filter(via_ecr == 0, arrival >= 0) %>%
  mutate(wait_time_call = case_when(
    !is.na(reneged) ~ 5,
    is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
    TRUE ~ call_taken - arrival # has been answered
  ))

baselineB_callwait_replication_summary <- baselineB_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(is.na(call_end) & is.na(reneged)),
            countfinished = sum(!is.na(call_end) | !is.na(reneged)))

# show top 10 results
head(baselineB_callwait_replication_summary, n = 10)
replication median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished
1 3.825883 3.183036 3.479728 0 5 1.464796 5 0 1726
2 3.744117 3.221769 3.258414 0 5 1.555424 5 0 1770
3 3.775021 3.116034 3.682359 0 5 1.206570 5 0 1772
4 3.952666 3.235613 3.446565 0 5 1.535425 5 0 1722
5 3.810451 3.201863 3.433046 0 5 1.446567 5 0 1736
6 3.550199 3.053321 3.688564 0 5 1.141632 5 0 1699
7 3.525437 3.005393 3.874786 0 5 1.005038 5 0 1652
8 3.808467 3.141373 3.616938 0 5 1.237069 5 0 1733
9 3.892609 3.168985 3.635974 0 5 1.225592 5 0 1750
10 3.623343 3.068928 3.714392 0 5 1.107349 5 0 1662
baselineB_callwait_summary <- baselineB_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "baseline") # to identify it later

baselineB_callwait_summary
median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished scenario
3.765179 3.147497 3.580269 0 5 1.334405 5 0 1722.8 baseline
baselineB_callwait_summary$mean_wait
## [1] 3.147497
baseline_callwait_summary$mean_wait
## [1] 3.147497

6.1 Part 3 Scenarios with option B

Trajectory must be included inside the lapply for replication

We can now analyse results

6.1.1 Count reneged

envs_scenarioB %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))
countreneged_mean
603.16

n=752

6.1.2 Count F2F

envs_scenarioB %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))
countF2F_mean
171.56

n=165

6.1.3 Call waits

scenarioB_callwait_results <- envs_scenarioB %>%
  get_mon_attributes() %>%
  pivot_wider(-c(time), names_from = key, values_from = value) %>%
  filter(via_ecr == 0, arrival >= 0) %>%
  mutate(wait_time_call = case_when(
    !is.na(reneged) ~ 5,
    is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
    TRUE ~ call_taken - arrival # has been answered
  ))

scenarioB_callwait_replication_summary <- scenarioB_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(is.na(call_end) & is.na(reneged)),
            countfinished = sum(!is.na(call_end) | !is.na(reneged)))

# show top 10 results
head(scenarioB_callwait_replication_summary, n = 10)
replication median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished
1 3.749845 3.209906 3.181770 0.0000000 5 1.464999 5 0 1666
2 3.932084 3.337168 3.053690 0.0000000 5 1.716950 5 0 1770
3 3.890046 3.282838 3.134920 0.0000000 5 1.540661 5 0 1789
4 3.764143 3.236900 3.134823 0.0000000 5 1.548852 5 0 1712
5 3.901525 3.279950 3.099232 0.0000000 5 1.651958 5 0 1723
6 4.060991 3.371249 3.040204 0.0111140 5 1.751085 5 0 1764
7 4.017161 3.352447 3.025200 0.0009612 5 1.774429 5 0 1763
8 3.904756 3.298098 3.105040 0.0000000 5 1.632506 5 0 1755
9 4.038012 3.353227 3.106330 0.0000000 5 1.716001 5 0 1764
10 3.684992 3.213686 3.210820 0.0000000 5 1.500257 5 0 1708
scenarioB_callwait_summary <- scenarioB_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "scenario") # to identify it later

scenarioB_callwait_summary
median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait countunfinished countfinished scenario
3.890618 3.289315 3.122619 0.001257 5 1.598594 5 0 1747.56 scenario
scenarioB_callwait_summary$mean_wait
## [1] 3.289315
scenario_callwait_summary$mean_wait
## [1] 3.289315